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Abstract 

Commensurate-incommensurate change on the one-dimensional 5 = 1 bilinear-biquadratic model 
(T~t(a) = Ylii^i • + a(Sj • Sj + i) 2 }) is examined. The gapped Haldane phase has two sub- 
phases (the commensurate Haldane subphase and the incommensurate Haldane subphase) and the 
commensurate- incommensurate change point (the Ameck-Kennedy-Lieb-Tasaki point, a = 1/3). 
There have been two different analytical predictions about the static structure factor in the neigh- 
borhood of this point. By using the S0rensen- Affleck prescription, these static structure factors are 
related to the Green functions, and also to the energy gap behaviors. Numerical calculations sup- 
port one of the predictions. Accordingly, the commensurate- incommensurate change is recognized 
as a motion of a pair of poles in the complex plane. 
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I. INTRODUCTION 



Commensurate-incommensurate (C-IC) transitions induced by frustration are important 
problems in many-body quantum spin systems. Among them, a C-IC change with an ex- 
citation gap is observed in one-dimensional (ID) quantum spin models.i^ This change is 
not a phase transition without an excitation gap. Whereas theories of C-IC transitions with 
no excitation gap (e.g. the Pokrovsky-Talapov transition 4 ) have been developed, those of 
C-IC changes for quantum systems have not been yet. In some classical systems, analytical 
approaches to the C-IC change have been discussed,— 1 ^ and then a random phase approx- 
imation approach has been succeeded phenomenologically. 7 However, on the one hand the 
ID frustrated Ising model for finite temperature cannot be mapped onto the ID quantum 
case, on the other hand the transfer matrix for the 2D Ising model on the triangular lattice 
is non-symmetric, thus its correspondence to the ID quantum case is not a simple problem. 
Therefore, an independent analytical research for the ID quantum C-IC change is needed. 

There are typical quantum models which show the C-IC change; the ID S — 1/2 next- 
nearest-neighbor (NNN) modeli and the ID S = 1 bilinear-biquadratic (BLBQ) model.— It 
is common between these models that the C-IC change occurs at the solvable point; the 
Majumdar- Ghosh point^ in the ID S — 1/2 NNN model and the Affleck-Kennedy-Lieb- 
Tasaki (AKLT) poinl^* 1 ^ in the ID S = 1 BLBQ model. These solvable points are called as 
the disordered point At the disordered point, the correlation length is the smallest and the 
ground state is described by the matrix product state . 8l9ll ° The correlation length and the 
incommensurate wave number are not differentiable at the disordered point, although they 
are continuous. The structure factor (the Fourier transform of the correlation function) varies 
from the 2D Ornstein-Zernicke type (the modified Bessel function) in the commensurate and 
incommensurate regions to the ID Ornstein-Zernicke type (the pure exponential function) 
at the disordered point.— 

Recently, in order to explain the C-IC change, some analytical studies have been proposed. 
Fath and Siito have suggested that the C-IC change occurs because of the existence of higher 
derivatives in an effective Lagrangian of the ID S = 1 BLBQ models On the other hand, 
one of us (KN) has discussed the static structure factor.— These studies show two candidates 
for the static structure factor, although they do not necessarily decide between them. 

By the way, S0rensen and Affleck (S-A) have studied two spin correlations and energy 
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FIG. 1: Ground state phase diagram of the S = 1 bilinear-biquadratic model. 

gaps between the triplet and singlet states under the open boundary condition by means 
of field theoretic approaches,— although they have not concerned with the C-IC change. 
Applying the S-A method to the C-IC problem, we can calculate some parameters included 
in the dynamical structure factor, i. e. the Green function. In our previous paper,— we have 
already found that the incommensurate wave number can be calculated by the energy gap 
of edge states. In this paper, we attempt to determine the Green function. After that, the 
relation between the singularities in the Green function and the incommensurability will be 
clear, and then we will obtain a unified view among commensurate and incommensurate 
behaviors. 

In this stage, we summarize some known properties of the ID S = 1 BLBQ model with 
the Hamiltonian; 

H(a) = ^{S, • S m + a(Si ■ S m ) 2 }. (1) 

i 

The ground state phase diagram of this model is shown in Fig. ^ This model is solvable 
at the AKLT point 9,10 a = = 1/3. The ground state is the Valence-Bond-Solid (VBS) 
state with the lowest excitation gap at the mode k — ir. One calls a phase, the ground 
state of which is a unique disordered ground state with a finite gap to the excited states, 
as the Haldane phase after Haldane's conjecture.^ This phase extends from the Takhtajan- 
Bubujian (TB) poin t 1 7i 18 i 19 a = -1 to the Uimin-Lai-Sutherland (ULS) poin t 2n i 21 i 22 a = 1. 
At the TB or ULS points, the BLBQ model is also solvable, and has the gapless ground 
state with the soft mode k = 0,tt or k = 0, ±27r/3,— respectively. For a < —1, there is 
the gapped dimerized (Dimer) phase,— 124125 whereas the gapless trimerized (Trimer) phase 
for a > 1.-2 5 - Between the AKLT point and the TB point, the lowest excitation has the 
wave number k = it, while the lowest excitations have the incommensurate wave number 
fcic, 2ir/3 < \kic\ < vr between the AKLT point and the ULS point." 2 ^ The wave numbers 
of the lowest excitations are different in these two regions, since the C-IC change occurs at 
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the AKLT point.^i The Haldane phase, therefore, has two subphases; the commensurate 
Haldane subphase for — 1 < a < 1/3 and the incommensurate Haldane subphase for 1/3 < 
a < 1. 

In addition, the VBS state becomes increasingly significant in connection with quantum 
entanglements . 27 ^ 8 ^ 29 The entanglements have a close relation to the matrix product state 
as well as the C-IC change. Therefore, it will be useful for an understanding of the quantum 
entanglements to investigate near the AKLT point, i.e. the C-IC change point. 

The organization of this paper is as follows. In the next section, we summarize essential 
points of the static structure factor concerning the C-IC change. The analyticity of the 
static structure factor explains that the change between branch points and a pole in the 
static structure factor corresponds to the C-IC change. In Sec. II 111 we discuss the relation 
between the edge states and the Green function on the basis of the S-A prescription. From 
the analysis of this section and Sec. HU we expect some behaviors of the energy gap of edge 
states. Before we study the energy gap of edge states numerically, we discuss the lattice 
effect in Sec. IIVI In Sec. El we confirm the gap behavior of edge states numerically, which 
is related to the Kennedy degeneracy— The last section gives a summary and a discussion. 

II. STATIC STRUCTURE FACTOR AND INCOMMENSURABILITY 

In our previous papers , 1 ? 1 ^ we have discussed the functional forms of the static structure 
factor concerning the C-IC change. Before studying the relation between edge states and 
the C-IC change, let us briefly summarize the essential points of the static structure factor 
about the C-IC change. 

A. Analyticity of the static structure factor 

From previous numerical results, especially in Ref. [llj], one can find the static structure 
factor in each region as follows: 

1. In the commensurate region (a <C ccd), 

S(q) oc J— 2 . (2) 
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(a) 



(b) • • 

FIG. 2: Typical branch cuts of f(z) = (z 2 - d)' 1 / 2 . (a) f(-z) = f(z). (b) f(-z) = -f(z). 

2. At the disordered point (a = «d), 

s(q) oc — ^ — 2- ( 3 ) 

gZ _|_ m l 

3. In the incommensurate region («>«])), 

S{q) oc 1 = + 1 ■ (4) 

V (? - <?ic) 2 + ™ 2 V (g + gic) 2 + 

However, one cannot connect these three expressions continuously. 

Considering an analytic continuation of real S(q) to the complex plane, we can discuss 
S(q) in the complex q plane. In terms of the singularity in the complex q plane, there are 
poles at the disordered point, in contrast to branch cuts in the other regions. 

In order to unify these three expressions, we reconsider the relation between a pole and a 
branch cut. Considering the next function, we can transform a pole into a branch cut, and 
vice versa, 

f(z) = (z 2 - d)- 1 / 2 , (5) 

where d is a real parameter. This function has two branch points. Typical branch cuts of 
f(z) are shown in Fig. |21 In the case of the branch cuts (a), which connect each of branch 
points to infinite distance, f(z) can be expanded in a Laurent series near z = 0, and then 
f(z) is found to be an even function f(—z) = f(z). On the other hand, in the case of the 
branch cut (b) which connects both of branch points, f(z) is an odd function f(—z) = —f(z) 
since f(z) can be expanded at infinite distance [see Appendix |X] in detail]. When d — 0, a 
simple pole appears in the case (b), whereas the branch cuts remain in the case (a). Thus 
we select the branch cut (b), and then deal with f(z) as an odd function. Then we find 
f(q — mi) satisfies 



f(q + mi) = f(q-mi), (6) 
f((-q) + rhi) = -f(q-rhi), (7) 



5 



where m is a real parameter. Note that q, q, and — q belong to the same Riemann sheet. 
The static structure factor must satisfy several physical requirements (PRI): 

1. [reality on the real axis] S(q) = S(q). 

2. [parity] S{q)=3{-q). 

3. [algebraic singularity] S(q) is an analytic function of a complex variable q except for 
several algebraic singular points. 

4. [analyticity on the real axis] Singular points and branch cuts must not cross the real 
axis. 

The above requirements represent properties of S(q) on a fixed a. In addition, 

5. [a-dependency near ajy] S(q) is an analytic function of a real parameter a in the 
neighborhood of «d- 

6. [property at qd] S(q) is described with two simple poles at the disordered point. 

On the basis of these requirements, we can obtain two possible candidates of the static 
structure factor near the disordered point: 

S'smg(g) =Af(q + mi)f(q-mi), (8) 

or 

S sing (q) = A—[f(q + ™) ^ /(? ^ (9) 

where real parameters A, fh, and d depend on a.— Figure H3 shows singularities of f(q — mi) 
when a) d < 0, b) d = 0, and c) d > 0. fh represents a distance between the real axis and the 
center of two branch points. Equations (jSJ) and (JSJ) tend to 1/q 2 at q — > oo limit. The pre- 
factor Ai) 1m in the difference type function © is determined so that S sing (q) = A/ (q 2 + rh 2 ) 
when d = 0. Equation (jHJ) is the same one which has firstly been proposed by Fath and 
Siito^ and the other (0) is discussed by KN.— We would like to clarify the behavior of the 
static structure factor S(q) by using another approach. In the following sections, we will 
investigate which is more appropriate structure factor either Eq. (jHJ or (jHJ). 
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FIG. 3: Singularities of f(q — rhi) when a) d < 0, b) d = 0, and c) d > 0. 
B. a dependency 

In addition, we can discuss how parameters m, d, and A depend on a. Since the correlation 
decays purely exponentially at the disordered point, we obtain (f = 0,m>0ata = «d- 
Generally, the requirement for the amplitude is A ^ since the correlation function becomes 
perfectly zero for A = 0. Near «d we then expect that d, in, and A can be expanded in a 
Taylor series: 

d = di(a — «d) + d 2 (a — «d) 2 + 0((a — «d) 3 ), (10) 
m = m + mi(a — «d) + — od) 2 ), (11) 

and 

A = A + A 1 (a-ct D ) + O((ct-ct D ) 2 ). (12) 

Besides, PRI-4 in Sec. HE] means that m > \/—d when d < 0. 

The incommensurate wavenumber gic = therefore, behaves as 

gic = V" - «D\Aii + ek(a - a D ), (13) 

in the incommensurate region, and gic = in the commensurate region. On the other 
hand, the correlation length £, which is related to the closest singular point to the real axis, 
behaves as 

C 1 ocm- V^d, (14) 

in the commensurate region, and 

r'ocm, (15) 

in the incommensurate region. 
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C. Numerical difficulties in dealing with the static structure factor 

The previous consideration results that the static structure factor should be Eq. (jHJ) 
or (J0J). To select one from two possibilities, we may calculate numerically the correlation 
function with the DMRG method. However, there are some difficulties in dealing with the 
static structure factor directly. We require 

1) to calculate a long range correlation near the disordered point since the incommensu- 
rate wave number is small, although the correlation length is short, 

2) to consider how to avoid edge effects, 
and also, 

3) to improve accuracy in calculating the correlation function, since the correlation func- 
tion is less accurate than the energy eigen values. 

Though it is indirect, there is another approach which uses the energy eigen values under 
the open boundary condition (OBC). This method has high accuracy even near the disor- 
dered point. In addition, small size data are important since we need to investigate poles 
far from the real axis. We only need to relate the energy eigen values to the static structure 
factor. 

In the next section, we will discuss the relation between the static structure factor and 
the energy eigen values under OBC, according to the S-A prescription.— 

III. EDGE STATES AND GREEN FUNCTION 

In this section, we discuss a Green function based on the S-A prescription^ [see Appendix 
IB] in detail]. 

A. Modified S-A prescription 

Now we consider a Green function G(q, k) which is the Fourier transform of G(x, r) in 
Euclidean space-time. The Green function determines various physical quantities, which 
contain a static structure factor S(q) and an energy gap of edge states. Between the Green 
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function G(q, k) and a frequency uo q (or an energy of a boson particle with a wave number 
q), the following relation is given in Appendix iBl 

G ^ K ) = T^TT2i ( 16 ) 

where k is a imaginary frequency The static structure factor is obtained by applying the 
Fourier transform of G(q, k), and then limiting as r — > 0; 

which recalls the original relation. One can show the correlation function from the Fourier 
transform of the static structure factor. 

(S x -S y ) = J (17b) 

Next, we shall examine the relation between the Green function and edge states. The 
edge states mean the triplet states and the singlet state under OBC. Among these states, 
there is a energy difference: 

AE S y = ^triplet — -^singlet- (18) 

The energy gap of edge states is connected with the Green function by the path integral 
method. The details are given in Appendix El Here, we only show the relation between the 
energy gap of edge states and the Green function: 

S eS = (-1) L A 2 S' X • S' L J dT X dr L ^G{q, ^^-H-ta-n), (19) 

where the left hand side of Eq. (fT9"j) means an effective action which is associated with 
an effective Hamiltonian, S e g = J drli^i an d A is an interaction parameter between the 
S = 1/2 edge spins S'^x and neighboring fields (p. The integral over n or tl provides a 
factor of 5(k). Thus we obtain^ 

A£ ST (L-1) = (-1) L A 2 / ^-G(q,K = 0)e^ L - 1 \ (20) 

J 27T 

Comparing Eq. (JT7J) with Eq. (|20|) . we see that Eq. (|2Ti|) is more manageable. The reason 
is that the integrand of Eq. (|2Uj) . i.e. the Green function, has poles, while that of Eq. (fTTj). 
i.e. the static structure factor, has branch points. 
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B. From static structure factor to Green function 



In Sec. HH we have discussed the static structure factor. We can apply a similar discussion 
to the Green function. Corresponding to the static structure factor, the Green function is 
permitted to have the following functional forms: 

K) = ^Tm^RFW^ ' (c/ Eq - m (21) 

or 

G sing (q, k) = G+ n Jq, k) + G; ing (q, «), (cf. Eq. ©) (22a) 



where 



G tinM «) = ~, TTTH^TT—^-, l\ - (22b) 



- u '- '' /,•-' (A/m) -'((r/ | mif-d)' 



G tmg(lj K ) (^sing(9; K )) nas the singularities only in the upper (lower) half g-plane. They 
satisfy 



G ing (?,«) = G sin g ( 23a ) 



G s1ng(-?> K ) = G s^ng(9, «), (23b) 

and 

G s1ng(9> ~ K ) = G^ng(?, «)■ (23c) 

In Appendix El we show that the static structure factor (Eq. (J0J)) is deduced from Eq. (J22|) . 

As well as the static structure factor, the Green function G(q, k) (both Eqs. (|21jl and 
(I22|) ) must satisfy the following physical requirements (PRII): 



1. [reality on the real axes] G(q, k) = G(q,K). 

2. [parity] G(q, k) = G(—q, k), G(q, k) = G(q, —k) 

3. [algebraic singularity] G(q, k) is an analytic function of complex variables q and k 
except for several poles. 

4. [analyticity on the real axes] Poles must not cross the real q and k axes. 
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5. [a-dependency near ao] G(q, k) is an analytic function of a real parameter a in the 
neighborhood of a^. 

However, Eq. (j2~Tj) is different from Eq. ()22|) when d = while the static structure factors 
(Eqs. (jHI) and 0) are the same (cf. PRI-6 in Sec. Ill A|) . Another difference is that in the 
limit q — > oo Eq. (J2*T|) behaves as g~ 4 while Eq. (|2*2*|) as g~ 2 . Hence, it is easier to distinguish 
Eqs. (pHj) and (J22J) clearer than Eqs. (JBJ) and © near the disordered point a = «d- 



C. energy gap of edge states 



On the basis of the above discussion with Eq. (J2DJ), the energy gap of edge states obtained 
from Eq. (J2~T|) is 

A 2 e -m(L-l) 



AE ST (L-1) = (-1) L X 2 



:sm(Vd(L-l) + (f)(m,d)), 



my/dy/m 2 + d 
l) L 2e-™ (L " 1) sm(\/d(L - 1) + 0(m, <£)), 



for d > (or a > «d), where 0(m, d) = tan l (yd/m), and also 
A£ ST (L - 1) = 



42 -m(L-l) 

1) L A 2 , , , , sinh(\/^d(Z, - 1) + (j)(m, d)), 



my —d\/m 2 + d 
= (-l) L Ie~™ (L - 1) sinh(V = d(i - 1) + 0(m, d)), 

for d < (or a < cud), where 0(m, d) = tanh _1 (v / — d/rh). 

On the other hand, the energy gap of edge states about Eq. (|2~2l is 

iA 2 e -m(L-l) 



ASgr(i-i; 

for a > «d, and also 

AE ST (L-V 



■1) L X 2 ' 



■Vd 



sm(\/d(L - 1)), 



(-1) L Ae- r ~ n{L - 1] sm(Vd(L - 1)), 



42 -m(L-l) 

-1) L A 2 „ sinh(v^d(L - 1)), 



-d 



;-l) L Ae-™ (L - 1) sinh(V = rf(i - 1)), 



(24a) 



(24b) 



(25a) 



(25b) 



for a < q;d- 

We will verify which is more appropriate between these two predictions (Eqs. (J24j) and 
fl2HJ)) by analyzing numerical data in Sec. |V| Note that Eq. (J23)l is apparently different 
from Eq. (|2l)j) when L — 1: Eq. (|23j) is always equal to zero, whereas Eq. (J2%Jl is nonzero. 
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IV. IMPLEMENTATION FOR LATTICE 



In this section, we consider an effect of the lattice structure. Equations (}2~2j) and are 
not equal to zero even when L is not an integer number, and therefore they are incompatible 
with the lattice structure. To include the lattice structure, we must require S(q) = S(q + 2n) 
and G(q, k) = G(q + 2tt, k). 

Now, we organize the new physical requirements (PRIII) for the static structure factor 
and the Green function, considering the lattice structure. PRIII from 1 to 5 are the same 
as PRI and PRII. We add the requirement of the periodicity to PRIII: 

6. [periodicity] S(q) = S(q + 2n) and G(q, k) = G(q, K + 2ir). 
From PRIII-6, we derive another physical requirement: 

7. [singularity in the Brillouin zone] There are only four singular points (poles or algebraic 
singularity) in the first Brillouin zone (— tt < < n). 

All the information needed for any problem can be determined in this zone. 

Then, we can construct some static structure factors and Green functions, satisfying these 
requirements, and we show them in Sees. IIV Al and IIV Bl 

A. Infinite sum version 

The easiest way is to consider the infinite sum of the translated singular parts. The static 
structure factor has the form as 

oo 

S(q)= ^2 S Bins (q + 2irj) + S ies (q), (26a) 

j=~oo 

and the Green function has 

oo 

G(q,K)= G aiDg (q + 2irj,K) + G Icg (q, k), (26b) 

j=-oo 

where both S s i ng (q + 2-af) and G sing (q + 2irj) represent shifted singular terms. S Teg (q) and 
G Teg (q, k) are regular functions in the whole q plane, such as 

oo 

s reg(q) = ^2 a i cos(Zg), (27) 

1=0 
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where a\ is a real number. In Eq. (|26[). the singular terms correspond to long-range behaviors 
in the real space, whereas the regular terms correspond to model-dependent short-range 
behaviors. 

Note that the infinite sum (J26|) for Eq. (J2J) or Eq. is divergent, whereas that for Eq. 
(JHJ) or Eq. (jUJ) is convergent. 

B. Sine wave version 

Alternatively, substituting a 2n or 47r-periodic function p(q) for q in S'(g) or G(q, k), we 
also obtain a 27r-periodic static structure factor S(p(q)) or a 27r-periodic Green function 
G(p(q), k), respectively. We impose some constraints on the periodic function p(q) to satisfy 



1. p(q) is a holomorphic function. 

2. p(q + 27r) = p(q) or p(g + 2n) = —p(q). 

3. =p(q) or p(-q) = -p(q). 

4. p(q) = p(q) or p(q) = —p(q). 

5. The inverse function p{q)~ l is a single- valued function in the first Brillouin zone — 7r < 
Kg < 7T. 

6. lim g _>oP(<?)/g = 1- 

The above requirements determine the distribution of zeros of p(q). From Weierstrass's 
theorem for infinite products 33 - and the above constraints, the function p(q) is determined 

as 



Replacing q in S(q) and G(q, k) by p(q), the static structure factor can be described as 



PRIII: 



p(q) =2 sin-. 



(28) 



S(q) = S sing (p(q)) + S rcg (q), 



(29a) 



and the Green function as 



G(q, «) = G sing (p(q), k) + G TCg (q, «). 



(29b) 
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C. Contour 



Corresponding to both the infinite sum version and the sine wave version, the contour C 
of the integral over q in Eqs. (I17|) and (|2(J|) is described in Fig. 0] Solid circles mean poles or 

o;;b 

■ ■ * Re q 

-271 -K n 2k 

FIG. 4: Contour C for the integral over q in Eqs. (jl7|) and (|2U|) . 

branches for the incommensurate case, and broken circles for the commensurate case. I, II, 
III, and IV represent the contours {q\(3lq : — 7r — > ir) PI (Qq = 0)}, {^(Kg = tt) fl (Sg : — > 
oo)}, {q\(JRq : 7r — > — 7r)n(Q : g = oo)}, and {q\(JRq = — 7i)n(^sq : oo — > 0)}, respectively. The 
contributions of II and IV cancel each other out because of the periodicity. The contribution 
of III can be ignored since S(q) and G(q, k) ^ q~ 2 as q — > oo. We, therefore, obtain that 
§ c = /j. As a result, the integral of the infinite sum Green function (Eq. (126bj) ) is equal to 
Eqs. (|2"4"j) and (|23|) . A similar discussion can be applied to the static structure factor. 

Note that the integral of the sine wave Green function (Eq. (j29bj)) is different from Eqs. 
(J21J) and ((23). We consider it in detail in Appendix iDl 



III 
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V. NUMERICAL ANALYSIS 



Our aim in this study is to decide between Eqs. (jHJ) and 0- In the previous section, 
each behavior of the energy gap of edge states has been expected from Eq. (jHJ) or 0. In 
this section, we, therefore, carry out the numerical calculation of the energy gap between 
the triplet and singlet states, and verify whether the results correspond to the predictions 
(Eqs. (|2~4"j) and (|25|) ) with the use of a nonlinear least-squares (NLLS) fitting program, which 
needs appropriate initial values. Applying the previous results,-^ we guess the initial values 
first. 



Although we have calculated the incommensurate wave number qte i n Ref- 15], its 



analytical reasoning was unclear. Also, we have not so far investigated m (the distance 
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between the real axis and the center of two singular points) and A (the amplitude of the 
energy gap). On the basis of the S-A prescription, we will calculate them in this section. 
We will also trace the singularities in the commensurate region. 



A. Surveys of edge states and incommensurability 



We treat the S = 1 BLBQ chain under OBC 



N 



(30a) 



hi — Si ■ Sj + i + a(Sj • Sj + i) 



(30b) 



where N is the number of the sub-Hamiltonian hi and a is the interaction constant of the 
biquadratic term. 

Note that N = L — 1, where L = {6, 7, • ■ ■ , 14} is the chain length. We can treat long 
chains (L > 14). However, their significant digit is smaller than that of short chains (L < 14) 
since their amplitudes of the energy gap are exponentially small near the AKLT point. Thus 
we treat up to L = 14. We exclude data smaller than L = 6 since the short-range behaviors 
are affected by mo del- dependent regular terms, i.e. S Teg (q) and G Teg (q,K) in Eqs. (j2EJ) and 



Since there are two edge S — 1/2 spin freedoms at the AKLT point (a = an), the 
following degeneracy occurs: 



(S = 1/2) ® (S = 1/2) = (S = 0) © (S = 1), 



(31) 



which reflects the Z 2 x Z 2 symmetr y 30 i 34 i 35 i 36 Therefore, the singlet-triplet energy gap (or 
the gap of edge states) 



AE ST (N, a) = triplet (iV, a) - E singlet (N, a), 



(32) 



is zero for all length spin chains at the AKLT point: 



AE ST (N,a D ) = 0. 



(33) 



Note that in the thermodynamic limit the triplet states and the singlet state also become 
degenerate in the whole Haldane phase (—1 < a < 1), and thus the amplitude of the gap of 
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edge states goes to zero as iV — > oo: 

lim AE ST (N,a) = 0. (34) 

To avoid confusion, we call the degeneracy at the AKLT point as the AKLT degeneracy. 

Numerical results of the gap of edge states are shown in Fig. For a^«D, the AKLT 
degeneracy breaks down. We see oscillating behaviors in the gap of edge states for a > a-Q. 
This phenomenon has been predicted from Eqs. (|24a|) and (J25a|) . Note that for « < «d the 
sign of the gap of edge states is different between even length chains and odd length chains 
because the parity of the bulk is different among these chains.— 

0.0001 
5e-05 



-5e-05 
-0.0001 

0.3 0.32 0.34 0.36 0.38 0.4 
a 

FIG. 5: Energy gaps of edge states AE$t = ^triplet — -^singlet as a function of a for various sizes. 




B. Initial guess 

1. incommensurate wave number 

Since the gap of edge states AE ST is a function of a and N, and it is oscillating in the 
incommensurate phase, we can find the relation between a and N, taking account of the 
condition Ai?sT = 0. Then we consider the n-th zero point of the singlet-triplet gap, 

AE ST (N,a n (N)) = 0. (35) 

If we adopt Eq. (|25a|) . i.e. AE$t(N) ~ sin(giciV), in the incommensurate region, we can 
relate the incommensurate wave number qic with iV as 

TTTl 

gicK(iV)) = — , (36) 
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where n = 1,2,3, ■■■ . We have already found that they are fitted by a universal curve like 
^/a — Q!d-— We show ql c /(a — ao) as a function of a — ojd in Fig. El These data fit well 
with the following equation: 



d(a) = qf c = d\(a — «d) + d 2 {a — «d) 2 + — 



«d) 3 ), 



where d x = 11.230 ± 0.010 and d 2 = -65.76 ± 0.83. 



ll 



Q 

a 
i 

a 



10 



0.005 0.01 0.015 0.02 0.025 
a-a n 



(37) 



FIG. 6: Dependence of the incommensurate wave number qic on a — «d- 

If we adopt Eq. (j24aj) . the corrections of 0{1/N) m.d\ and d 2 should be found. Since their 
corrections are smaller than 2 %, we see that our guess adopting Eq. ()25a|) . AEst(N) ~ 
sin (VdN), is more reliable than Eq. ()24a)l . We expect that the gap of edge states behaves 
in the incommensurate and commensurate regions as Eq. (J25aj) and Eq. (j25b|) . respectively. 
This result means that the static structure factor is Eq. Q. 

In addition, the number of zero points n (except the AKLT point) and the system size 
N(= {1, ■ • • , 13}) are correlative. It is easy to find the following relation: 

n = N (mod 3). (38) 

We see that the relation between the max number of zero points n max and iV as 

(39) 



7m max (iV) 7i 
N < 3' 



We can confirm this relation up to = 13 (n max (13) = 4). We expect that 7m max (A r )/A^ 
has the limit n/3 as N tends to oo, and therefore the position of the max-n-th zero point 
a n max goes to the ULS point as iV — > oo. 
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2. amplitude and center of coupling poles 



In the previous sub-subsection we have found that the Green function corresponds to the 
difference type of the static structure factor (Eq. (JSJ)). Next, we examine the parameter A 
and fh. 

In the incommensurate region a > «d, we expect that the gap of edge states has the 
following form: 



AE ST (N) 



\N+1 



Ae~ mN sin(gicJV). 



(40) 



The incommensurate wave number q\c{ot) near the AKLT point is obtained in the previous 
sub-subsection. Using these values and considering the following equation, we can determine 
A and m: 



log 



AE ST (N) 



log I A I — mN. 



(41) 



sin(giciV) 

Figure [7| shows log \AEst/ sin(giciV)| when a = 0.3492 behaves linearly as a function of N. 
We have just confirmed our prediction for the incommensurate region. 




10 11 12 13 



FIG. 7: Finite size results for log | AE ST / sin(gi C iV)| when a = 0.3492. 

Similar consideration can be applied to the commensurate region. In the region, the gap 
of edge states should be 

AE ST (N) = (-l) JV+1 le-^sinh(g c ^), (42) 

where qc = y/—d = a/— (a — ar>)\fdi — d 2 (a — «d)- Here, the commensurate wave num- 
ber qc is indirectly determined by using d\ and d 2 , which are obtained from Eq. (pTTj) . 
Considering the following equation, we can obtain A and rh in the commensurate region: 

AE ST (N) 



log 



sinh(g c iV) 



log \A\ — mN. 



(43) 
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C. Nonlinear least-squares fitting 



In the previous subsection, we have adopted the commensurate wave number qc indirectly 
determined by using the parameters of the incommensurate wave number in Eq. (J37j) . 
although the relation between qic and qc is somewhat unclear. Actually, it seems that the 
region, where Eq. ()37|) is permitted, may be narrower in the commensurate side than in the 
incommensurate side. In order to check the above mentioned results from another viewpoint, 
we employ a NLLS fitting method. Using this method, we can determine above parameters 
A, in, and d directly. Since the method requires appropriate initial values, we must have 
determined them in the previous subsection. 

Taking into account the fact that the amplitude of the energy gap is exponentially small 
near the AKLT point, we use the following weighted values to perform the NLLS fitting 
program: 

y N = (-l) N+1 AE ST (N)w N , (44) 

where wn = exp(m'iV) is a weight, and ih' = ih + 5 is a value estimated from ih at the 
nearest a. Correctly, what we determine by the NLLS fitting method is not ih but 5. 
The NLLS fitting method requires the minimization of the squared residuals, 

Q = '52\(VN-M*)) i , ( 45 ) 

where x = (A, ih, d) and /jv(x) is a fitting function of x. From the minimum value of Q, we 
obtain parameters (A,ih,d). 

In the case of a - «d = 0.02, for example, we show the data of the energy gap and 
the fitting function /(x) = A exp(— ihN) sm(^/dN) where A = —0.421, ih = 0.991, and 
d = 0.202 in Fig. |H1 Also, the case of a — «d = 0.1 is shown in Fig. |H1 The parameters of 
fitting function /(x) are A = -0.758, ih = 0.638, and d = 0.810. 

1. fitting with Eq. (|25|) 

Figure El summarizes the fitting results with Eq. (J25|) in the incommensurate region. 
The obtained parameters, A, ih, and d for < \a — «d| < 0.05, in which region Q is less 
than 1.0 x 10 -8 , are shown in Figs. ^2 E2 and 1131 respectively. Near the AKLT point, 
they converge very well, and behave continuously with a. 
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FIG. 8: Fitting results with Eq. (J25J) when a - a D = 0.02. 
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FIG. 9: Fitting results with Eq. (J2HJ) when a - «d = 0.1. 

We see that these parameters are smooth between the commensurate and incommensurate 
regions. In Fig. ^2 ^ is highly linear. In Figs. ^] and each parameter m, d varies 
linearly in the incommensurate region, whereas there are a broad maximum and a broad 
minimum, respectively, at a — «d — —0.02 in the commensurate region. The range where 
d, in, and A can be expanded in terms of a — «d is narrower in the commensurate region 
than in the incommensurate region. When a — «d is less than -0.02, there should be a 



different mechanism from what we have expected in Sec. Ill Bl since PRIII-5 is not satisfied 
in the region. 

Now, we estimate the correlation length £ from the obtained data. Usually, the correlation 
length is related to an inverse of a distance between the closest singular point and the real 
axis. In the incommensurate region, £ = m" 1 , while £ = (m — V - i n the commensurate 
region. These results are shown in Fig. This behavior is consistent with the previous 
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FIG. 10: Nonlinear least-squares fitting results with Eq. (|25(l . 
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FIG. 11: Nonlinear least-squares fitting results for A. 

numerical result.— 
2. fitting with Eq. 

We also attempt to apply the nonlinear least-squares fitting program to Eq. (J2H). The 
obtained parameters (A, m, d) are shown in Fig. ED We see that they behave as discontin- 
uous pieces about a. In addition, the region where d oc a — old is very narrow. These facts 
mean that the supposed functions (Eq. (jHJ) and Eq. ()24|)) are not correct. Of course, the 
residual Q is larger than the one shown in the previous sub-subsection. 
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FIG. 12: Nonlinear least-squares fitting results for in. 
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FIG. 13: Nonlinear least-squares fitting results for d. 
3. about sine wave version 

We have so far assumed the infinite sum version (Eq. (J26j) ) as a lattice effect. Now we 
consider the case of the sine wave version (Eq. ()29j) ). The energy gap of edge states (in the 
incommensurate region) is modified as 

AE ST {L - 1) = (-1)^1^(0(^-1) Bin(3(C)(£ - 1)), (cf. Eq. flliD) (46) 

where ( = — 21og(— iz + y/l — z 2 ) and z = (mi + \fd)/2 [see Appendix ID], The obtained 
parameters with the NLLS fitting program are shown in Fig. In this figure, the param- 
eters behave continuously except for some discontinuous points near a — «d = 0.025 and 
0.065. Comparing Figs. E3andE3 we think that the results of the infinite sum version (Eq. 
()26|0 as a lattice effect is more reasonable than that of the sine wave version (Eq. (j2Hl))) 
although we have not found a conclusive evidence to support it yet. 
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FIG. 14: Correlation length for < \a - «d| < 0.05. 
1.5 
1 

^0.5 
S 

< 
-0.5 
-1 

0.02 0.04 0.06 0.08 0.1 

a-a D 

FIG. 15: Nonlinear least-squares fitting results with Eq. I)24JI . 

VI. SUMMARY AND DISCUSSION 

In this study, we have examined the S = 1 BLBQ model near the AKLT point. Analyzing 
the energy gap of edge states on the basis of the S-A prescription,— we have shown that 
our numerical results support Eq. ()25|) . i.e. Eq. © which is one of the predictions for the 
static structure factor concerned with the C-IC change. The energy gap of edge states is 
more manageable than the correlation function because the singularities are different among 
them, and thus our results are clearer than the previous one. We have also obtained the 
incommensurate wave number, the amplitude and the correlation length. These results 
are consistent with the previous result.— Our incommensurate wave number qic is different 
from the original incommensurate wave number kic in Sec. U] The difference is caused of 
our notation; the prefactor (— 1) 7V+1 of the gap Ai?sT is left apparently. Two different wave 
numbers can be related as kic — ft ± Qic- 
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FIG. 16: Nonlinear least-squares fitting results with Eq. (|4(ijl . 

We should mention here that Eq. (|23j) is not only numerically supported, but also it 
has a physically favorable feature. From Eq. (J23j) . one can see A.Est(0) = for L = 1, 
i.e. N = which means no sub-Hamiltonian case in Eq. (}30|) . Although this property is 
not necessary since overall G(q, k) consists of singular and regular terms as Eqs. (j26b|) and 
(|29bj) . the property A£?st(0) = seems quite natural physically 

The amplitude A has been found to be proportional to \/a — «d- This result implies that 
A 2 oc a — «d because of Eq. (|23|) . In Appendix |B] and Ref. [l^, we have only assumed that 
the interaction A is some real constant. However, our results suggest that A is some complex 
constant. Thus we have to modify the assumption for A. Note that A is equal to zero at the 
disordered point, corresponding to the VBS picture. 

Originally in the S-A prescription, the singlet-triplet energy gap AE^t depends on the 
Green function, which is assumed to have a simple pole in the upper half-plane and in the 
lower half-plane. However, our results suggest that two poles should be concealed in the 
upper or lower half-plane. In general, one of these poles is far from the real axis, and therefore 
the ordinary field theoretic approach, like the nonlinear a model,-^ appears to succeed in 
describing the Haldane phase. Indeed, if we explain the whole Haldane phase including the 
C-IC change, we must consider the four singular points. Near the AKLT point, a four-pole 
structure becomes explicit in the Green function, and then the incommensurability occurs 
in the incommensurate Haldane subphase. A prelude to the incommensurability arises even 
in the commensurate region. We have found that positions of poles (singularities) included 
in the Green function are represented in terms of (rh, d). 
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We have left some future tasks; the effective Lagrangian (maybe two components) and 
the dispersion curve for the Green function Eq. (}2*2*|) (cf. those for Eq. (}2"Tj) have been 
obtained in Ref. numerical verification of the static structure factor and the dynamical 
structure factor. Although we treat only the ID S — 1 BLBQ model in this paper, we have 
obtained similar results about the ID S = 1/2 NNN model.— However, we need to modify 
the discussion about the Green function since a quasiparticle has a magnon-like behavior in 
5=1 models, whereas a spinon-like behavior in 5* = 1/2 models. 
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APPENDIX A: DOUBLE- VALUED FUNCTION f(z) 

In this appendix El we examine some properties of f(z) (Eq. (jSJ)) introduced in Sec. 



1. Choice of branch cuts and related property 

The function 

f(z) = (z 2 - d)' 1 ' 2 = (z + sfd)-V 2 (z - Vd)^ 2 (Al) 

is a double- valued function with two branch points at z = —\fd and z = \fd. We can 
freely choose branch cuts of f(z) although the parity of the selected branch cut should be 
compatible with that of f{z). Typical branch cuts are shown in Fig. El (a) both of the 
branch points are connected, and (b) each of them are connected to infinite distance. These 
different branch cuts bring different parities to f(z). 
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FIG. 17: Two typical types of branch cut of f(z) = (z 2 — d) 1 ^ 2 . The points A and B show 
z = —\fd and y/d, respectively. 

We consider the case (a) first. We can carry out the Laurent expansion of f(z) around 
z = oo when \z\ > -\/\d\'- 

/2 1 ~ (2»-l)M /rf\° 

(z - d> = z^-^ar\?) ' (A2) 

It is an odd function with a zero point of order one at infinity. 

About the case (b), we can expand f(z) around z = when \z\ < \f\d\, and then we 
obtain an even function: 

l'-^=(^iyg)" (A3) 

Alternatively, we can explain their different parities by a graphical way. Let z + Vd = re ia 
and z — Vd = pe %13 . Then 

f( Z ) = r -l/2 /0 -l/2 e -i(a+/3)/2_ (M) 

In the case (a) 

f(- Z ) = r -l/2p-l/2 e -j(7r+/3+7r+a)/2 

= -/(*), (A5) 

and in the case (b) 

f(- Z ) = r -l/2 p -l/2 e -i(/3-7r+7r+a)/2 

= f(z). (A6) 
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2. First and second sheets 



The Riemann surface of f(z) consists of two Riemann sheets. Here, we consider a relation 
between the first and second Riemann sheets (z\- and 22-plane, respectively), although we 
focus on the case that both of branch points are connected by a branch cut. As shown in 
Fig. UHl let (1 and (2 be a point on the z\- and Z2-plane, respectively, although these two 
points have the identical coordinate. A similar discussion with the previous subsection can 
be applied to the case of (1 — ■> (2. Then we find 



/(Ca 



-/(& 



(A7) 




Zj Z 2 



o o 



FIG. 18: Two Riemann sheets of f(z) = (z 2 - d)- 1 / 2 . Ci is a point on the first sheet (zi-plane), 
and £2 is on the second sheet (z2-plane). 



APPENDIX B: FIELD THEORETIC APPROACH FOR EDGE STATES 

In this appendix[B]we reproduce the S0rensen and Affleck prescription. 14 They start from 
the nonlinear a (NLu) models Since an effective field model is not clear in our case, it is 
not possible to apply this model as it is near the AKLT point. However, we can develop a 
similar discussion if we assume a Green function G(q, k). The Green function is determined 
from discussions in Sees. ITT1 and ITTT1 It describes some massive free boson fields <f>(x,r): 

4>(x,t) = [ r ^ <1 / (a(g)e i9X -^* + a 1 ^ (q) e - iqx+i ^ t ) , (Bl) 

J y2'K^2u q 

where a(g) is a bose operator which satisfies [a(g), a^(g')] = 8(q — q'), and u q is obtained 
from the Green function. Vacuum expectation values among two different boson fields are 
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calculated as the following: 



J iqx-iw q t 

(4>(x,t)- 0(0,0))= / — . (B2) 



27T 2u q 



1. Static structure factor 



Static structure factor S(q) is defined as the Fourier transform of a equal-time correlation 
function: 

(<t>(x,0) >(0,0)) = | ^(g)e^. (B3) 
Therefore, we can relate the static structure factor with uj q : 



Si,) = (B4) 



2. Green function 



Green function is defined as the time-ordered expectation value: 

iG(x, t) = T(<t>(x, t) ■ 0(0, 0)). (B5) 

Using the Wick rotation 

t = —IT, UJ = —IK, (B6) 
where ut = —kt, and the step function 

«M = S J da l^' (B7) 

we then find 

f C iqX+iKT iC(q k) — f 1 ( e iqx+i{-iu q +a)T j c -i(-iu> q +a)T-iqx\ 

J (2tt) 2 ' J 2(27r) 2 iau q 



dndq eiqx+iKT i 

2 ,,-2 I ,.,2 ' 



(2tt) 2 k 2 + ujI 

where k is a imaginary frequency. The Green function G(q, k) associates with u: q as 



(B8) 



G{q,K) = -±— 2 . (B9) 
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3. Perturbation theory 



Using the Green function G(q, k), we can describe a free part of the action. 

S&) = \ J ntfe*"*", (BIO) 

where (f> is the Fourier transform of <fi. 

Next, we take the edge effects into consideration. The open boundaries have the effect of 
leaving a S = 1/2 degree of freedom at each end of chain. The edge spins will interact with 
the rest of the system. To consider this effect, we assume the following interaction: 

Hi = A[0(1) • S\ + {-l) 1 - 1 ^) ■ S' L ), (Bll) 

where A is weak coupling constant. S'i and S'^ are two S = 1/2 excitations known to exist 
at the end of the open chai n 30 ' 34 ! 35 ' 36 . The sign in front of the second term comes from the 
reason that we consider the boson field <fi with the wave number tt. 

Carrying out the ordinary Gaussian integral, we can obtain an effective action 
Seff(S i, S'l). 

J J)(j) e ~ S i ( t')+I dTdxJ(x,T)-ct>(x,T) _ (J e ~ S cS^ (B12) 

where J(x, r) = X[S'i8(x — X\) + (—1) l ~ 1 S'l5(x — Xl)]- Then we find 

S eS = (-l) L \ 2 S\ ■ S' L I dndnj0pG(q, K ) e ^ L -^^\ (B13) 

The constant C in Eq. (jB12|) contains the divergent self-energy that comes from terms 
with both arguments included in the Green function on the same source world-line. These 
correspond to virtual 4> particles that are emitted and absorbed by the same source. We are 
not interested in these, but only in the variation in the vacuum energy as a function of the 
separation of the sources. 

In this appendix, we have not consider an imaginary time dependency of S'i 5 l since such 
a dependency has so far been unclear. 

APPENDIX C: TRANSFORMATION FROM GREEN FUNCTION TO STATIC 
STRUCTURE FACTOR 

In this appendix O we will show that the static structure factor (Eq. @) is constructed 
from the Green function (Eq. (|22p). 
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We consider the following integral: 

f e iKT dn ine iTy ^ 

J (K - y/z) (K + y/z) y/z 

where r > 0, z = re td , and 9 = Argz (0 < 9 < 2ir). The right hand side of Eq. (jC^lj) is a 
double- valued function and has a branch point at z — 0. 

Now we consider w = y/z. In general, w corresponds to W\ = y/re tQ l 2 in the upper half 
w-plane when < argz < 2tc, while w corresponds to W2 = — y/re l6 l 2 in the lower half 
if-plane when 2n < arg z < 4tt. Thus Eq. (jQTj) is rewritten as 



(C2) 



e iKT dK \me lTW1 lw\ (0 < argz < 2n) 



[k y/z)(K + y/z) \ i7[e ir W 2/ W2 (2vr < arg 2 < 4vr). 
Using Eq. (|C2J) . we can show that 



/dx An 
—e^(G + (q, K ) + GL(g, «)) = — (/(g + mi) - /(? - mi)). (C3) 

APPENDIX D: INTEGRATION OF GREEN FUNCTION ABOUT SINE WAVE 
VERSION 

Substituting p(g) = 2sin(g/2) for q in Eq. (|22|). we obtain the energy gap behavior of 
the edge states; 

AE ST (L - 1) = (-1) L A 2 4 / ^ L ^{G+(p(q), K ) + G-(p(g), «)} (Dl) 

mr J 2ir 

Using the formula sin -1 z = ilog(—iz + \/l — z 2 ), we can integrate the right hand side of 
Eq. (|D1|) over q. After the integration, we find 



AE ST (L - 1) = (-l)^ e W-i) sin(3(C)(L - 1)), (D2) 
where ( = —2 log(— iz + \/I — z 2 ) and z = (mi + y/d)/2. 
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